function fQ = check_fQ(grid,base,f,Q,Pi)
% fQ =
%
% Compute the integral of f*Q which, with proper boundary conditions,
% should be zero.

 fQ = zeros(grid.d,1);
 for ie=1:grid.ne
   for l=1:base.m
     Qv(:,l) = base.o(:,:,l)*Q(:,ie);
   end
   Qv = grid.e(ie).b*Qv / grid.e(ie).det_b;

   fQ = fQ + ...
       grid.e(ie).det_b * sum( ...
            base.wg                                    ...
        .* (f(grid.e(ie).b*base.xig+grid.e(ie).x0)-Pi) ...
        .* Qv                                          ...
                             ,2);
   
 end

%
% for id=1:grid.d
% odi = permute(base.o(id,:,:),[2,3,1]);
%   for ie=1:grid.ne
%     odi_p = grid.e(ie).b * permute(base.o(id,:,:),[2,3,1]) /grid.e(ie).det_b
%
%     fQ(id) = fQ(id) + ...
%       grid.e(ie).det_b * sum( ...
%            base.wg                                    ...
%        .* (f(grid.e(ie).b*base.xig+grid.e(ie).x0)-Pi) ...
%        .* (Q(:,ie)'*odi)   ...
%                             );
%   end
% end

return

